home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / DMFSS.FOR < prev    next >
Text File  |  1985-11-29  |  7KB  |  246 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE DMFSS
  5. C
  6. C        PURPOSE
  7. C           GIVEN A SYMMETRIC POSITIVE SEMI DEFINITE MATRIX ,DMFSS WILL
  8. C           (1) DETERMINE THE RANK AND LINEARLY INDEPENDENT ROWS AND
  9. C               COLUMNS
  10. C           (2) FACTOR A SYMMETRIC SUBMATRIX OF MAXIMAL RANK
  11. C           (3) EXPRESS NONBASIC ROWS IN TERMS OF BASIC ONES,
  12. C               EXPRESS NONBASIC COLUMNS IN TERMS OF BASIC ONES
  13. C               EXPRESS BASIC VARIABLES IN TERMS OF FREE ONES
  14. C           SUBROUTINE DMFSS MAY BE USED AS A PREPARATORY STEP FOR THE
  15. C           CALCULATION OF THE LEAST SQUARES SOLUTION OF MINIMAL
  16. C           LENGTH OF A SYSTEM OF LINEAR EQUATIONS WITH SYMMETRIC
  17. C           POSITIVE SEMI-DEFINITE COEFFICIENT MATRIX
  18. C
  19. C        USAGE
  20. C           CALL DMFSS(A,N,EPS,IRANK,TRAC)
  21. C
  22. C        DESCRIPTION OF PARAMETERS
  23. C           A     - UPPER TRIANGULAR PART OF GIVEN SYMMETRIC SEMI-
  24. C                   DEFINITE MATRIX STORED COLUMNWISE IN COMPRESSED FORM
  25. C                   ON RETURN A CONTAINS THE MATRIX T AND, IF IRANK IS
  26. C                   LESS THAN N, THE MATRICES U AND TU
  27. C                   A MUST BE OF DOUBLE PRECISION
  28. C           N     - DIMENSION OF GIVEN MATRIX A
  29. C           EPS   - TESTVALUE FOR ZERO AFFECTED BY ROUND-OFF NOISE
  30. C           IRANK - RESULTANT VARIABLE, CONTAINING THE RANK OF GIVEN
  31. C                   MATRIX A IF A IS SEMI-DEFINITE
  32. C                   IRANK = 0 MEANS A HAS NO POSITIVE DIAGONAL ELEMENT
  33. C                             AND/OR EPS IS NOT ABSOLUTELY LESS THAN ONE
  34. C                   IRANK =-1 MEANS DIMENSION N IS NOT POSITIVE
  35. C                   IRANK =-2 MEANS COMPLETE FAILURE, POSSIBLY DUE TO
  36. C                             INADEQUATE RELATIVE TOLERANCE EPS
  37. C           TRAC  - VECTOR OF DIMENSION N CONTAINING THE
  38. C                   SOURCE INDEX OF THE I-TH PIVOT ROW IN ITS I-TH
  39. C                   LOCATION, THIS MEANS THAT TRAC CONTAINS THE
  40. C                   PRODUCT REPRESENTATION OF THE PERMUTATION WHICH
  41. C                   IS APPLIED TO ROWS AND COLUMNS OF A IN TERMS OF
  42. C                   TRANSPOSITIONS
  43. C                   TRAC MUST BE OF DOUBLE PRECISION
  44. C
  45. C        REMARKS
  46. C           EPS MUST BE ABSOLUTELY LESS THAN ONE. A SENSIBLE VALUE IS
  47. C           SOMEWHERE IN BETWEEN 10**(-4) AND 10**(-6)
  48. C           THE ABSOLUTE VALUE OF INPUT PARAMETER EPS IS USED AS
  49. C           RELATIVE TOLERANCE.
  50. C           IN ORDER TO PRESERVE SYMMETRY ONLY PIVOTING ALONG THE
  51. C           DIAGONAL IS BUILT IN.
  52. C           ALL PIVOTELEMENTS MUST BE GREATER THAN THE ABSOLUTE VALUE
  53. C           OF EPS TIMES ORIGINAL DIAGONAL ELEMENT
  54. C           OTHERWISE THEY ARE TREATED AS IF THEY WERE ZERO
  55. C           MATRIX A REMAINS UNCHANGED IF THE RESULTANT VALUE IRANK
  56. C           EQUALS ZERO
  57. C
  58. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  59. C           NONE
  60. C
  61. C        METHOD
  62. C           THE SQUARE ROOT METHOD WITH DIAGONAL PIVOTING IS USED FOR
  63. C           CALCULATION OF THE RIGHT HAND TRIANGULAR FACTOR.
  64. C           IN CASE OF AN ONLY SEMI-DEFINITE MATRIX THE SUBROUTINE
  65. C           RETURNS THE IRANK X IRANK UPPER TRIANGULAR FACTOR T OF A
  66. C           SUBMATRIX OF MAXIMAL RANK, THE IRANK X (N-IRANK) MATRIX U
  67. C           AND THE (N-IRANK) X (N-IRANK) UPPER TRIANGULAR TU SUCH
  68. C           THAT TRANSPOSE(TU)*TU=I+TRANSPOSE(U)*U
  69. C
  70. C     ..................................................................
  71. C
  72.       SUBROUTINE DMFSS(A,N,EPS,IRANK,TRAC)
  73. C
  74. C
  75. C        DIMENSIONED DUMMY VARIABLES
  76.       DIMENSION A(1),TRAC(1)
  77.       DOUBLE PRECISION SUM,A,TRAC,PIV,HOLD
  78. C
  79. C        TEST OF SPECIFIED DIMENSION
  80.       IF(N)36,36,1
  81. C
  82. C        INITIALIZE TRIANGULAR FACTORIZATION
  83.     1 IRANK=0
  84.       ISUB=0
  85.       KPIV=0
  86.       J=0
  87.       PIV=0.D0
  88. C
  89. C        SEARCH FIRST PIVOT ELEMENT
  90.       DO 3 K=1,N
  91.       J=J+K
  92.       TRAC(K)=A(J)
  93.       IF(A(J)-PIV)3,3,2
  94.     2 PIV=A(J)
  95.       KSUB=J
  96.       KPIV=K
  97.     3 CONTINUE
  98. C
  99. C        START LOOP OVER ALL ROWS OF A
  100.       DO 32 I=1,N
  101.       ISUB=ISUB+I
  102.       IM1=I-1
  103.     4 KMI=KPIV-I
  104.       IF(KMI)35,9,5
  105. C
  106. C        PERFORM PARTIAL COLUMN INTERCHANGE
  107.     5 JI=KSUB-KMI
  108.       IDC=JI-ISUB
  109.       JJ=ISUB-IM1
  110.       DO 6 K=JJ,ISUB
  111.       KK=K+IDC
  112.       HOLD=A(K)
  113.       A(K)=A(KK)
  114.     6 A(KK)=HOLD
  115. C
  116. C        PERFORM PARTIAL ROW INTERCHANGE
  117.       KK=KSUB
  118.       DO 7 K=KPIV,N
  119.       II=KK-KMI
  120.       HOLD=A(KK)
  121.       A(KK)=A(II)
  122.       A(II)=HOLD
  123.     7 KK=KK+K
  124. C
  125. C        PERFORM REMAINING INTERCHANGE
  126.       JJ=KPIV-1
  127.       II=ISUB
  128.       DO 8 K=I,JJ
  129.       HOLD=A(II)
  130.       A(II)=A(JI)
  131.       A(JI)=HOLD
  132.       II=II+K
  133.     8 JI=JI+1
  134.     9 IF(IRANK)22,10,10
  135. C
  136. C        RECORD INTERCHANGE IN TRANSPOSITION VECTOR
  137.    10 TRAC(KPIV)=TRAC(I)
  138.       TRAC(I)=KPIV
  139. C
  140. C        MODIFY CURRENT PIVOT ROW
  141.       KK=IM1-IRANK
  142.       KMI=ISUB-KK
  143.       PIV=0.D0
  144.       IDC=IRANK+1
  145.       JI=ISUB-1
  146.       JK=KMI
  147.       JJ=ISUB-I
  148.       DO 19 K=I,N
  149.       SUM=0.D0
  150. C
  151. C        BUILD UP SCALAR PRODUCT IF NECESSARY
  152.       IF(KK)13,13,11
  153.    11 DO 12 J=KMI,JI
  154.       SUM=SUM-A(J)*A(JK)
  155.    12 JK=JK+1
  156.    13 JJ=JJ+K
  157.       IF(K-I)14,14,16
  158.    14 SUM=A(ISUB)+SUM
  159. C
  160. C        TEST RADICAND FOR LOSS OF SIGNIFICANCE
  161.       IF(SUM-DABS(A(ISUB)*DBLE(EPS)))20,20,15
  162.    15 A(ISUB)=DSQRT(SUM)
  163.       KPIV=I+1
  164.       GOTO 19
  165.    16 SUM=(A(JK)+SUM)/A(ISUB)
  166.       A(JK)=SUM
  167. C
  168. C        SEARCH FOR NEXT PIVOT ROW
  169.       IF(A(JJ))19,19,17
  170.    17 TRAC(K)=TRAC(K)-SUM*SUM
  171.       HOLD=TRAC(K)/A(JJ)
  172.       IF(PIV-HOLD)18,19,19
  173.    18 PIV=HOLD
  174.       KPIV=K
  175.       KSUB=JJ
  176.    19 JK=JJ+IDC
  177.       GOTO 32
  178. C
  179. C        CALCULATE MATRIX OF DEPENDENCIES U
  180.    20 IF(IRANK)21,21,37
  181.    21 IRANK=-1
  182.       GOTO 4
  183.    22 IRANK=IM1
  184.       II=ISUB-IRANK
  185.       JI=II
  186.       DO 26 K=1,IRANK
  187.       JI=JI-1
  188.       JK=ISUB-1
  189.       JJ=K-1
  190.       DO 26 J=I,N
  191.       IDC=IRANK
  192.       SUM=0.D0
  193.       KMI=JI
  194.       KK=JK
  195.       IF(JJ)25,25,23
  196.    23 DO 24 L=1,JJ
  197.       IDC=IDC-1
  198.       SUM=SUM-A(KMI)*A(KK)
  199.       KMI=KMI-IDC
  200.    24 KK=KK-1
  201.    25 A(KK)=(SUM+A(KK))/A(KMI)
  202.    26 JK=JK+J
  203. C
  204. C        CALCULATE I+TRANSPOSE(U)*U
  205.       JJ=ISUB-I
  206.       PIV=0.D0
  207.       KK=ISUB-1
  208.       DO 31 K=I,N
  209.       JJ=JJ+K
  210.       IDC=0
  211.       DO 28 J=K,N
  212.       SUM=0.D0
  213.       KMI=JJ+IDC
  214.       DO 27 L=II,KK
  215.       JK=L+IDC
  216.    27 SUM=SUM+A(L)*A(JK)
  217.       A(KMI)=SUM
  218.    28 IDC=IDC+J
  219.       A(JJ)=A(JJ)+1.D0
  220.       TRAC(K)=A(JJ)
  221. C
  222. C        SEARCH NEXT DIAGONAL ELEMENT
  223.       IF(PIV-A(JJ))29,30,30
  224.    29 KPIV=K
  225.       KSUB=JJ
  226.       PIV=A(JJ)
  227.    30 II=II+K
  228.       KK=KK+K
  229.    31 CONTINUE
  230.       GOTO 4
  231.    32 CONTINUE
  232.    33 IF(IRANK)35,34,35
  233.    34 IRANK=N
  234.    35 RETURN
  235. C
  236. C        ERROR RETURNS
  237. C
  238. C        RETURN IN CASE OF ILLEGAL DIMENSION
  239.    36 IRANK=-1
  240.       RETURN
  241. C
  242. C        INSTABLE FACTORIZATION OF I+TRANSPOSE(U)*U
  243.    37 IRANK=-2
  244.       RETURN
  245.       END
  246.